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Theory of quasiparticle spectra and the de Haas- van Alphen (dHvA) oscillation in type-II super- 
conductors are developed based on the Bogoliubov-de Gennes equations for vortex-lattice states. As 
the pair potential grows through the superconducting transition, each degenerate Landau level in 
the normal state splits into quasiparticle bands in the magnetic Brillouin zone. This brings Landau- 
level broadening, which in turn leads to the extra dHvA oscillation damping in the vortex state. We 
perform extensive numerical calculations for three-dimensional systems with various gap structures. 
It is thereby shown that (i) this Landau-level broadening is directly connected with the average gap 
at H = along each Fermi-surface orbit perpendicular to the field H; (ii) the extra dHvA oscillation 
attenuation is caused by the broadening around each extremal orbit. These results imply that the 
dHvA experiment can be a unique probe to detect band- and/or angle-dependent gap amplitudes. 
We derive an analytic expression for the extra damping based on the second-order perturbation with 
respect to the pair potential for the Luttinger-Ward thermodynamic potential. This formula repro- 
duces all our numerical results excellently, and is used to estimate band-specific gap amplitudes from 
available data on NbSe2, NbaSn, and YNi2B2C. The obtained value for YNi2B2C is fairly different 
from the one through a specific- heat measurement, indicating presence of gap anisotropy in this 
m aterial. C programs to solve the two-dimensional Bogoliu bov-de Gennes equations are available 
at tittp: //phys . sci .hokudai . ac . jp/~kita/index-e . html. 



YNi 2 B 2 
URu 2 Si 2 1111 




I. INTRODUCTION 

The de Haas-van Alphen (dHvA) experiment on nor- 
mal metals has been .a_iinique and powerful tool to probe 
their Fermi surfaces.En3 The main purpose of this paper 
is to establish theoretically that it can even be used to 
detect detailed gap structures of type-II superconductors. 

Back in 1976, Graebner and Robbins discovered the 
dHvA oscillation in 2H-NbSe 2 persisting down through 
the superconducting upper critical field i? c2 J3 It was af- 
ter. 15 years later when Onuki et al. first reconfirmed 
ito Since then, however, a considerable number of ma- 
terials have been found to display the dHvA oscilla- 
tion in the VjOEjtex state™ They include: A15 supercon- 
ductors VaSjata NbaSnjD a borocarbide superconductor 
heavy- fermion superconductors CeRu 2 ,E£I 
UPd 2 Al 3 J13 CeCoIn 5 ,£il and.an organic su- 
perconductor k-(BEDT-TTF) 2 Cu(NCS) 2 ;E£| see Refs. p| 
and [ll] for a recent review. Basic features of the oscil- 
lation are be summarized as follows: (i) the dHvA fre- 
quencies remain unchanged through the transition; (ii) 
the oscillation amplitude experience an extra attenua- 
tion; (iii) the cyclotron mass does not change except for 
strongly correlated heavy fermion systems. 

It is somewhat surprising that the dHvA oscillation 
is observable even in superconductors without a well- 
defined Fermi surface. Many theories have been pre- 
sented to expl^iiLJj^pipislcua^ and the ex- 
tra damping,ymaBmeBBBeEl which may 
be classified into three categories. 

The first approach applies a Bohr-Sommcrfcld semi- 
classicali— iquantization to either the Braadt-Pesch- 
Tewordty (BPT) Green's function near H c2 }B the elec- 
tron number N at H = ,E_£l or Dyson's equation at 
H = 0,c3 for obtaining the oscillatory behavior of the 



magnetization. As may be seen by this diversity of the 
applications, however, there is no unique semiclassical 
quantization scheme for quasiparticles superconductors; 
thus, the validity of the procedure is not clear. This cat- 
egory includes Maki's theory,E3 whictt^was later repro- 
duced by Wasserman and Springfordcil by treating the 
BTP self-energy as the extra broadening factor in the 
normal-state thermodynamic potential and then follow- 
ing Dingle's procedure.E3 However, the BTP self-energy 
itself is obtained within the quasiclassical approximation 
without the Landau-level structure so that this approxi- 
mation may also be questionable. 

The second approach relies on some approximate an- 
alytic solutions for the Bogoliubov-de Gennes (BdG) 
equations or the equivalent Gor'kov ecura tions, such as 
the coherent-potential ax>proximation,E§E3 the diagonal- 
pairing approximation^ or the Ginzburg-Landau (GL) 
expansion far the free energy with respect to the order 
parameter.t3n3o However, quantitative estimations of 
those approximations are yet to be carried out. It should 
also be noted that the GL expansion necessarily inte- 
grates out the quasiparticle degrees of freedom, so that 
the physical origin of the extra oscillation damping may 
be obscured in the GL approach. 



The third approach solves thepJMG equations 
merically without approximations.E3c3 Norman et a/.K 
thereby extracted an analytic formula for the dHvA os- 
cillation damping through a fitting to their numerical 
dataO However, it is a two-dimensional calculation for 
the isotropic s-wave pairing where the number of Lan- 
dau levels below the Fermi level is N-p ~ 10 at iJ c2 . 
As may be realized from the appearance of the quan- 
tized Hall effect, two-dimensional systems in high mag- 
netic fields may be qualitatively different from three- 
dimensional systems. Thus, the obtained formula may 
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not be appropriate for describing real three-dimensional 
materials with Np ^> 1. On the other hand, another 
calculation by— Miller and Gyorffy for a two-dimensional 
lattice modeEJ corresponds to the low-field limit near 
i? c itij and may not be suitable to explain the experi- 
ments. Moreover, calculations for lattice models have a 
flaw that they cannot yield continuous magnetic oscilla- 
tion due to the commensurability condition between the 
underlying lattice and the vortex lattice. 

Notice finally that most of the above theories consider 
only the isotropic s-wave pairing. Especially, no numeri- 
cal studies have been performed yet for anisotropic pair- 
ings or three-dimensional systems. 

With these observations, we will perform both nu- 
merical and analytic calculations for three-dimensional 
BdG equations with various gap structures. They can be 
solved efficiently by the Landau-level-expansion method 
(LLX), which was formulated-for vortex-lattice states of 
arbitrary pairing symmetrycfl and used successfully to 
compare low-enetgy quasiparticle spectra between s- and 
d-wave pairings £3 We will thereby clarify how the dis- 
crete Landau levels experience quantitative changes as 
the pair potential grows below H C 2- Another purpose 
is to find out the connection between the gap anisotropy 
and the extra dHvA amplitude attenuation. Terashima et 
aln reported a dHvA experiment on YM2B2C where an 
oscillation is observed to persist down to a field ^Q.2H C 2- 
On the other hand, a specific-heat experiment at H =dl 
shows a power-law behavior ocT 3 at low temperatures,E£l 
indicating existence-of gap anisotropy in this material. 
Indeed, Izawa et alxJl recently reported presence of four 
point nodes in the gap-based on a thermal-conductivity 
measurement. MiyaketJ argued that point or line nodes 
along the extremal orbit may weaken the damping, and 
proposed to use the dHvA effect as a probe for gap 
anisotropy. We examine this possibility in full details 
and present a quantitative theory on the issue. 

This paper is organized as follows. Section II provides 
a formulation to solve the BdG equations for vortex- 
lattice states. Section III presents calculated quasiparti- 
cle spectra and the dHvA oscillation for two-dimensional 
systems to clarify their basic features as well as the origin 
of the extra dHvA oscillation damping in the vortex state. 
In Sec. IV, it is demonstrated that the gap anisotropy 
at H = can be detectable via the dHvA oscillation in 
the vortex state based on both numerical and analytic 
calculations for three-dimensional systems with various 
gap structures. Section [v| presents estimations of energy 
gap for NbSe2, NbsSn, and YNi2B2C using the analytic 
formula obtained in Appendix Section VI concludes 
the paper with a brief summary. In Appendix |a], we 
derive a convenient expression for the thermodynamic 
potential. Appendix B summarizes the expressions of 
basis functions and overlap integrals used in the numer- 
ical calculations. In Appendix C, we derive an analytic 
expression for the extra dHvA oscillation damping in the 
vortex state. A brief report of the contents was already 
presented in Ref. pq. 



II. FORMULATION 

A. Bogoliubov-de Gennes equations 

Throughout the paper we will rely on the mean-field 
BdG equations, which obtain the quasiparticle wavefunc- 
tions u s and v* with a positive eigenvalue E s > by 



dr 2 



H(r 1 ,r 2 ) A(n,r 2 ) 
At(n,r 2 ) -2f(n,r 2 ) 



u s (r 2 ) 
- v J( r 2) 



E s 



u s (ri) 
-v s *( ri ) 



(1) 



Here A is the pair potential and H denotes the normal- 
state Hamiltonian in the magnetic field; both are 2x2 
matrices to describe the spin degrees of freedom. The 
symbol t denotes Hermitian conjugate in both the coordi- 
nate and spin variables as [AJ(ri, r 2 )] CTlC r 2 = A* a2t7i (r 2 , r±) 
with cTj =|, J,. With this definition, we can see immedi- 
ately that the matrix in Eq. ([!]) is Hermitian. 
We adopt the free-particle Hamiltonian for H- 



2£(ri,r 2 ) = <S(ri-r 2 ) 



[-zW 2 + fA(r 2 )]' 
2m e 




where m e , —e (e > 0), c, and £f are the electron mass, 
the electron charge, the light velocity, and the chemi- 
cal potential, respectively. We will not consider the spin 
paramagnetism throughout. We also neglect the spatial 
variation of the magnetic field as appropriate for the rel- 
evant high-K materials. Then, the vector potential A can 
be expressed using the symmetric gauge as 



A(r) 



--Bzxr , 
2 



(3) 



where B denotes the average flux density, and we have 
chosen the field along — z for convenience. 

The pair potential in turn is given with respect to the 
quasiparticle wavefunctions as 



A(ri,r 2 ) = y(ri-r 2 )$(ri,r 2 ), 



(4) 



where V denotes the interaction and $ is the order pa- 
rameter defined by 



<£(ri,r 2 ) 



E 



[u s (n)vj(r 2 ) - v s (ri)uj(r 2 )] 



1 , E s 
x - tanh — — — 

2 2k B T 



(•5) 



with T the temperature and T denoting the transpose. 

It is shown in Appendix ^ that the thermodynamic 
potential corresponding to Eqs. (|l])-(^|) is given by 

Q = -fc B T^ln(l+e-^/ fcBT ) ~J2 E * [ Ws(r)\ 2 dr 

s s ^ 

- X - I / , TrA t (r 1 ,r 2 )$(r2,r 1 )dr 1 dr 2 , (6) 
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where Tr denotes taking trace over spin variables. This 
expression will be useful to obtain an analytic expression 
for the extra dHvA amplitude attenuation in the vortex 
state. The magnetization is then calculated by 



M 



g(n/v) 

OB 



(7) 



where V is the volume of the system. 



B. Vortex lattices and basic vectors 

Solving the above equations for general non-uniform 
systems is a formidable task. For lattice states, however, 
it can be reduced into a numerically tractable problem 
using the symmetry that they are periodic with a single 
flux quantum </>o = hc/2e per unit cell. We hence define 
a pair of basic vectors by 



'ai = (ai x ,aiy,0) 
a 2 = (0,02,0) 



(ai xa 2 ) 



90 
B 



ttI% , 



(8) 



Here ipNka is a quasiparticle basis function with N de- 
noting the Landau level, k defined by Eq. (|io|), and a 
(=1,2) signifying signifying two-fold degeneracy of every 

orbital state. On the other hand, ip NcCl and W Ntm are 
basis functions for the center-of-mass and relative coor- 
dinates, respectively, with iV c and N r denoting the cor- 
responding Landau levels, q defined by Eq. (pT|), and m 
an eigenvalue for the relative orbital angular momentum 
operator l z . The quantities p z and L are, respectively, 
the wavenumber and the system length along the z di- 
rection parallel to the magnetic field; we adopt a notation 
of using p as a wavevector in zero field to distinguish it 
from the two-dimensional magnetic Bloch vector k per- 
pendicular to the field. As noted in Ref. [39[ an arbitrary 
single q suffices to describe the conventional Abrikosov 
lattices due to the broken translational symmetry of the 
vortex lattice. Then the first expansion of Eq. ( |l2| ) tells 
us that, by choosing q = 0, we get an complete analogy 
with the uniform system in that the pairing occurs be- 
tween (k,p z ) and (— k, — p z ). Finally, the two expansion 

coefficients A^|^ 2 and A^,™' 1 are connected by 



where Ib = y/hc/eB is the magnetic length. The basic 
vectors of the corresponding reciprocal lattice are then 
defined by 



b 1 = 2(a 2 xz)//| 
b 2 = 2(zxai)//| 



(9) 



We now introduce .magnetic Bloch vectors for quasipar- 
ticle eigenstates byE.il 



bi 



Aft Aft 



- Aft 1 Aft z 
and those for the center-of-mass coordinate by 



Ml 



bi 



Aft Aft 



Aft . Aft 



(10) 



(ii) 



where Aft is an even integer with Aff denoting the number 
of flux quanta in the system. Notice that q covers an area 
four times as large as that of k. 



C. Landau-level-expansion method (LLX) 

It has been shownfEl that the pair potential of the con- 
ventional Abrikosov lattice can be expanded in two ways 
with respect to (i"i,r 2 ) and (R,r) = ( ri l" r2 , ri— r 2 ) as 



A(ri,r 2 ) 



p ipz(zi— z 2 ) 



E E ^NiNi ^JVika(ri) ^JV 2 q-ka(r2) 
ka NiN 2 



= 7=E E (-i«f*(R)^! m w : 



Nc N I mp z 



N c N r m 



x( _ 1) ^ r ) j 



(13a) 



A^ = jjr £ (A c A r |7V 1 7V 2 )^(m + A r |2k- q > 

f N t N 2 k 



(13b) 



where (NiNt\N c N T ) and (2k— q\m+N T ) are elements of 
unitary matrices for the basis change, i.e. overlap inte- 
grals. Their explicit expressions together with those for 

tpNka, and ^wlm are S iven in Appendix ||. 

A great advantage of using Eq. ( |l2| ) is that it enables 
us to transform Eq. (Q) into a numerically tractable prob- 
lem, as mentioned before. Indeed, expanding the quasi- 
particle wavefunctions as 



u(r)= ^ Us (N)4> Nka {r)—=, (14a) 

Nkap z V 

v(r) = J2 v ^) ^ q -ka(r) , (14b) 

Nkap z * 



Eq. ([l]) is reduced to a separate matrix eigenvalue prob- 
lem for each kap z , and the eigenstate is labelled by 
s — v\<iap z a with v and a denoting the quasiparticle band 
and its spin, respectively. Explicitly, Eq. (|l|) becomes 



E 



oy(P*) \{^Pz) 
LLn 1 N 2 —N ± N 2 

A( k P^)t _n-/(P") 
i±N 1 N 2 J^N-i_N 2 



(12) where A^ p ^ is given by Eq. (13a), and HniN2 ^ s diag 



" u s {N 2 )~ 




= E S 


. _-v:(^ 2 )_ 



u s (Ax) 



(15) 



/(p.) 
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onal as 



ear in the pair potential! 



<iy(Pz) _ r 



(Ni + \)hw B + 



2m n 



I, (16) 



with WB = eB / Mc the cyclotron frequency. 

The self-consistency equation (||) are also simplified 
greatly. Let us define 

V pp , = J ^(r) C - i (P- p ') r d 3 r 

oo t 
£=0 m=-l 

where limfp) = ®em(dy>)-i^= exp(imcp p ) is the spherical 

harmonic O We also expand both A and $ in terms of 
the center-of-mass and relative coordinates as the last 
line of Eq. (|l2|). Then Eq. (Q) is transformed into an 
equation for the expansion coefficients of each (N c , m) as 



1 



2ttEL 



J2 Vi(p,p')e em (6 p )B im {6 p ,)&£ s 



(N c m) 



(18) 



with p = V n JI'b+pI and 6 p = tan" 1 ( VN ^ 1b ) . 

Let us further assume that a single £ is relevant in Eq. 
( |l7| ) and take the corresponding Ve(p,p') in a separable 
form as 



V t (p,i/) = ViWt(£)Wt(?), 



(19) 



(JV c m) 



2 ^ 

JViJV 2 



(NcNrlNiNz) 



tanh & +tanh If 



2T 



27 



^(-lr^^lAT.+niV.-^AX^ 



6+6 

{N c +n m+n) 



(21) 



Equation ( p0| ) with Eq. (|T]) determines the mean-field 
H C 2(T), i.e. T C (H). If we use the asymptotic expres- 
sion (B6) for the overlap integral (N c N r \NiN2) and re- 
place the sum over Ni by the integral over x — (Ni — 
N2)/ y/2(N c + N T ), we reproduce the smooth quasiclassi- 
cal _ff^ asl (T) obtained, for example, for the s-wave pair- 
ing by Helfand and Werthamer.Cj 

Finally, Eq. (^) for two-dimensional systems can be 
transformed similarly. It is also obtained from Eqs. (|l7|)- 
( pp| ) by replacing Ve(p,p') — * (p,p r ), extending the 
summation over m in Eq. ( |l7| ) from —00 to 00, and finally 
restricting the summation over I and p z only to £ = and 
Pz = 0, respectively. 



III. TWO-DIMENSIONAL CALCULATIONS 

We first consider a couple of two-dimensional mod- 
els and perform fully self-consistent calculations. Our 
purposes in this section are summarized as follows: (i) 
to clarify the essential features of the results by self- 
consistent calculations; (ii) to see whether point nodes 
in the .gap really enhances the dHvA signals as Miyake 
claims .113 



where Wt (£) is some cut-off function with respect to £ = 
h 2 p 2 /2m r — e-p satisfying W^(0) = 1. Then can rewrite Eq. 
( |l8| ) as 



A 



A {Ncm) w e (ti)®im(e P ) 



-N T p 

with £ = h 2 (N T /1%+pl) /2m e - e F and 



^(JV c m) 



Vt 



2-kIIL 



Y^Wt(z')e im (o P ')k 



(N c m) 

p' ) 2tNy M 



(20a) 



(20b) 



A. Models 

The one-particle Hamiltonian (|J) yields an isotropic 
Fermi surface specified by a unit vector _p = 
(cos tp p , sin if p ). As for the pairing interaction (|l7|), we 
adopt the following models: 



^pp' 



V W(QW{?) 

V2w®(pi-i$w(e)(p»-p*) 



, (22) 



Thus, we only need a self-consistent solution for a set of 

" ™\t,B)} through Eqs. (|Tf 

) using Eq. 



and 

It has been shownE3EJ that retaining a few iV c 's, e.g. 
iV c = 0, 6, 12 for the hexagonal lattice, is sufficient to de- 
scribe the vortex lattices of H>0.05H C 2- Thus, the orig- 
inal problem of obtaining self-consistency for A(ri,r2) 
at all space points is now reduced to the one for a few 

expansion coefficients {A*' c (T, B)}. This situation is 
analogous to the zero-field case where a single parameter 
A (T) specifies the pair potential. 

The linearized self-consistency equation is obtained by 



substituting into Eq. (20b) the expression of J> 



lin- 



where 



W(0 = exp 



2 V huJD 



(23) 



is a smooth cut-off function with a cut-off frequency. 
The second model of Eq. (^||) is beyond the original 
isotropic interaction (|l7|), but it is convenient for the 
above-mentioned purposes. In zero field, the two inter- 
actions yield the s- and d,, 2^2 -wave gaps as 



- p ~\A Q W(0(pl~P 2 y ) l ci2 



s-wave 



d x 2_ y 2-wave 



(24) 
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respectively. The corresponding A^™" 1 of Eq. ( 53 ) for 
B || z is given by 



r(JV c m) 



(25a) 



with t; = h 2 N T /2m c lg — eF and 
' Vb 



aw = 



j^E^K')*?'" 



^2 



*W>.2) ,^(iV ,-2) . (25b) 
AT' ' Ni 



Here we have adopted a normalization for A^-* different 
from Eq. (poh so that acquires a direct correspon- 

dence to the maximum gap Ao in Eq. (pi|); the factor i 
in the second case stems from cos 2(p p in Eq. (|2~5 
We have chosen Vg in Eq. (p2|) as 



g t = -N(0)V e = 0.5, 



(26) 



where A^(0) = m e /2irH 2 is the density of states per spin 
at the Fermi level. Another important parameter is the 
zero-temperature coherence length defined by 



(27) 



with vp the Fermi velocity. We have adopted pf£,o =5 for 
our calculations. The above two quantities fix our models 
completely; the cut-off Huif> in Eq. ( p3| ) and T c (B = 0) are 
calculated using the gap equation. 

It should be noted that choosing pf£,o also determines 
the following quantities: (i) the ratio hw HV iB*i / k-QT Cl 
where fej^quasi is the zero-temperature cyclotron en- 
ergy at the quasiclassical upper critical field -ff^ 2 uasl ; (ii) 
the number Nf of the Landau levels below the Fermi 
level at ff^ asl . Indeed, using the usual cut-off model 
W(^) = 9(hiOF> — |£ P |) with 8 the step function, we obtain 
the following results for the s-wave pairing: 



/kJ^quasi/fceTc = 6.28/pf£o > 

N F = e F /hw H ^ = 0.140(p F £o) 2 . 



(28) 



To reproduce the experimental situation tno H q_^i / k^,T c — 
1^3 within the present model, we should go into the 
quantum limit pf£o = 6^2, but we then have Nf = 5 ~ 
1. In real materials, however, Huj ffqU a B i / 'k^T c and pf£o 
are apparently independent parameters due to effects not 
covered by the free-particle model such as the energy 
band structure. Indeed, pf£o is of the order of 30 (NbSe2) 
or even larger, whereas Hu> H quasi /k^T c = 1 ~ 3; see Table |l| 

below. Also, Nf ^S> 1 for those materials. The failures to 
describe these situations are among the main difficulties 
of the free-particle model of Eq. (p|) . 

Motivated by these observations, we also perform an- 
other calculations with much more Landau levels below 



the Fermi level. This is achieved by including the ef- 
fect of the band dispersion. Specifically, we apply the 
Onsager-Lifshiz (OL) quantization scheme to H of Eq. 
([!]), i.e. the procedure which has been very successful £p* 
describing the dHvA oscillations in the normal stateBEl 
Given the density of states per spin N(e) and the average 
flux density B, the Nth Landau level En (N — 0, 1, 2, ■ • • ) 
is determined by 







K 


if— 2 / 







N(e') de' . 



(29) 



WefixW^ 



of Eq. ( |15| ) in this way assuming it is diag- 
onal. In contrast, we use the same expression for A.^^ 
of Eq. (15) as the free-particle case. Finally, we choose 



£ = £N t /2-£F in Eq. 
free-particle model.l 



25 ) based on the consideration of the 
With these prescriptions together 
with the transformation (p^|), the coupled equations ( |l5| ) 
and (|2^) are defined unambiguously. We adopt the model 
density of states: 



N(s) 



2ttH 2 



1 



aT 



r 2 



(30) 



and choose the numerical constants (a, T) = (2.1, 2.7) and 
(5.0, 1.0) for the s- and d-wave models of Eq. (p4|), respec- 
tively. We also use Eq. (^6| ) for the pairing interaction 
and fix fiion = 0.5ef in Eq. (]23|). We thereby obtain 



k B T c at T = and N F ~ 30 at H = H 



T^quasi 



which describe the experimental situation much better 
than the free-particle model. 



B. Numerical procedures 



Coupled equations (15) and (|25j) are solved iteratively 
with the help of the transformation (|l^) to obtain self- 
consistent {A^'j-'s and the corresponding quasiparticle 
eigenstates. The hexagonal (square) vortex lattice is as- 
sumed for the s-wave (d g 2_ y 2-wave) .model, as expected 
theoretically in high magnetiG-EeldsSa and observed re- 
cently in Lai.83Sro.i7Cu04+,5.Ej It should be noted, how- 
ever, that the precise lattice structure is not important 
for the theory of the dHvA oscillation in superconductors. 
We set q= ^(bi+b2) in the relevant equations so that_a. 
core of the pair potential is located at the origin R = 0.l2I 
An advantage of this choice is that the corresponding 
quasiparticle energies have the rotational symmetry of 
the hexagonal (square) lattice around k = 0; other choices 
would shift the rotation axis from the orig in. We then 
perform calculations of Eqs. (|l5|) and (|25| ) for a set of 
discrete k's defined by Eq. (|10[), where A/f is chosen as 
a multiple of 12 to include all the high- symmetry points 
r, M, and K (T, X, and M) of the hexagonal (square) 
lattice. Three different values A/f = 12, 24, 36 are used to 
see the size dependence, and it has been checked that the 
results do not differ for the three cases. The hexagonal 
(square) symmetry enables us to restrict the summation 
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FIG. 1: The upper critical field H r p. as a function of T for 
(a) s-wave and (b) d-wave of Eq. (p4j). Here Pf£o = 5, and 
H C 2 is normalized by the quasiclassical upper critical field 
# c q 2 uasi (T = 0). 




r m k 



o 
S3 




r x m 



over k into approximately 1/12 (1/8) area of the Brillouin 
zone. Thus, the calculations can be reduced greatly with 
due care on the degeneracy of high-symmetry points. Fi- 
nally, the obtained eigenvalues and eigenstates are sub- 
stituted into Eq. (^) to calculate the magnetization by 
Eq. (0). All the calculations are performed at T = 0.1T C . 




j3 



0.8 1.0 




0.8 1.0 

B/HT%0~) 




0.8 1.0 

B/ffr"(o) 



FIG. 2: The expansion coefficients A'^- 1 in Eq. (^5|) as a 
function of B for the s-wave (first column) and the d-wave 
(second column) with T = 0.1T C and pf£o = 5. The dotted 
lines in the first row signify the square-root behavior expected 
from the quasiclassical theory. 
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FIG. 3: Quasiparticle dispersion in the magnetic Brillouin 
zone for the s-wave hexagonal lattice (first column) and the d- 
wave square lattice (second column) . HerepF^o = 5, T = 0.1T C , 
and B/if c q 2 uasi (0) is equal to 0.8, 0.5, and 0.1 from top to 
bottom, respectively 



C. Results 

The above self-consistent procedure is known to give 
rise to oscillatory singular behaviors in both H C 2 and 
range where the dHvA oscillation 
persfetsBBfia Fi gure [l] displays H c2 (T) calculated self- 
consistently for the s- and <i-wave models; it is normalized 
by the quasiclassical upper critical field _£Z^ 2 uasl (T = 0). 
An oscillatory behavior sets in around T < 0.2T C , and 
H C 2 deviates substantially from the smooth Helfand- 
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FIG. 4: Oscillatory part of magnetization M osc for (a) the s- 
wave and (b) the d-wave, over 0.8 < asi (0)/B < 2.0 (1.2 > 
B/i^ 2 uasi (0) >0.5) with T = 0.1T C and p F £o = 5. The dotted 
lines are the curves of the corresponding normal state. 
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Werthamer behavioru3 predicted by the quasiclassical 
theory. The number of the Landau levels below the 
Fermi level is Ap ~ 10 around -ff^ asl (0), which is con- 
siderably smaller than those for the real materials. Fig- 
ure | shows A (N ^ as a function of B at T = 0.1T C for 
the s-wave hexagonal lattice (first column) and the d- 
wave square lattice (second column); they are real and 
finite only for A c = JLrL2 , • • • (0, 4, 8, • • • ) for the hexag- 
onal (square) lattice,EJe2l as already mentioned. We ob- 
serve that A( Nc ''s are also singular, and the dominant 
A(°) component cannot be described by the square-root 
behavior near H^ 2 uasl (0) expected from the quasiclassi- 
cal theory. However, those singular behaviors disappear 
gradually as B decreases. 

Figure ^| displays the quasiparticle energies in the mag- 
netic Brillouin zone for the s-wave hexagonal lattice (first 
column) and the d-wave square lattice (second column) 
at B/iJ^ asl (0) = 0.8, 0.5, and 0.1. At 5/ff c q 2 uasi (0) = 0.8, 
we already observe large dispersion for E< 2Aq where the 
pair potential is effective. In contrast, the flat Landau- 
level structure remains for E > 2Ao where the pair po- 
tential vanishes in the present cut-off model of Eq. (23). 
Thus, the dispersion is caused clearly by the scatter- 
ing from the growing pair potential, and as will be dis- 
cussed below, it is the origin of the extra dHvA oscillation 
damping in the vortex state. We also notice that, for 
B/H^2 ml (0) <; 0.5, almost no qualitative difference can 
be seen between the s- and d-wave cases. At a lower field 
oi B/H^ 2 uasl (0) — 0.1, however, a marked difference grows 
around E < Ao. The s-wave energy bands of E < 0.7Ao 
are flat and occur in pairs with the level spacing of the o*-i 
der of Ap/ep- As already pointed out by Norman et aZ.,E3 
these corresponds to the bound core states of an isolated 
vortex with little tunneling probability between adjacent 
cores. In contrast, the d-wave bands in the same re- 
gion are densely packed with large dispersion, indicating 
the extended nature of the corresponding quasiparticle 
wavefunctions. From this comparison, we conclude that 
no bound states exist for the d-wave model even in the 
zero-field limit of an isolated vortex, in agreement with 
the result of Franz and Tesanovic.caThis difference in the 
low-energy dispersion at low fields was already reported 
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FIG. 5: The expansion coefficients A^ Nc ^ in Eq. ( p5| ) as a 
function of B for the s-wave (first column) and the d-wave 
(second column). Here T = 0.1T C , and the non-quadratic dis- 
persion given by Eq. dio| ) is used. The dotted lines in the 
first row signify the square-root behavior expected from the 
quasiclassical theory. 



in Ref. p5| 

Figure 14 shows oscillatory part of the magnetization 
Af osc calculated numerically by Eq. (0), where curves of 
the corresponding normal state are also plotted for com- 
parison. The damping starts from above 7?^ asl (0) where 
A(°) is already finite as in Fig. ||, and develops rapidly 
as grows in decreasing B. Thus, the mean-field the- 
ory predicts that the dHvA oscillation comes together 
with the oscillatory singular behaviors in H C 2 and A( Wc '. 
Combined with the energy dispersion given in Fig. [|, we 
are now able to attribute the origin of the extra damp- 
ing unambiguously to the Landau-level broadening due 
to the pair potential. The oscillations are rather irregular 
in both the s-wave and d-wave cases, in accordance with 
the singular behaviors of A^- 1 in Fig. |[ We also see no 
qualitative difference between the two cases. 

However, the free-particle model has several inappro- 
priate points as discussed already around Eq. (|2^). For 
example, the number of Landau levels below the Fermi 
level Af is necessarily Af ~ 10 at id^ a81 for pf£o = 5, 
which is much smaller than the values of the materi- 
als displaying the dHvA oscillation. Hence the above 
numerical results may not be sufficient to say anything 
quantitative about the dHvA attenuation or the differ- 
ences between the s- and d-wave cases. We have thus 
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FIG. 6: Oscillatory part of magnetization M oac for (a) the 
s-wave and (b) the d-wave, over 0.8 < i?^ lasl (0)/-B < 1.7, 
i.e. 1.25 > B/H^ si (0) > 0.59. Here T = 0.1T C , and a non- 
quadratic dispersion given by Eq. (^o|) is used. The dotted 
lines are the curves of the corresponding normal state. 



FIG. 7: Oscillatory part of magnetization A/ osc for (a) the 
s-wave and (b) the d-wave. The difference from Fig. M lies 
in the use of quasiclassical A^'s given by the dotted lines in 
Fig. | 



performed another calculations for the model described 
around Eqs. (|29| ) and (|3Cj ) where huj H q Ua si /ksT c ~ 1 and 

7V F ~30at 5 = # c q 2 uasi (0). 

Figure || shows the field-dependence of the expansion 
coefficients As Nc > calculated self-consistently for the s- 
wave hexagonal lattice (first column) and the d-wave 
square lattice (second column). Singular oscillatory be- 
haviors are manifest in both cases as in the case of the 
quadratic dispersion, which originate from the singular 
density of states of Landau levels.E3 For example, the 
dominant A^ ) component have a nonzero value from 
above H^^ 1 and deviates substantially from the qua- 
siclassical square- root behavior (dotted lines). 

Figure || displays the corresponding oscillatory part of 
the magnetization M osc calculated numerically by Eq. 
(0), where normal-state results (dotted lines) are also 
plotted for comparison. The main features are summa- 
rized as follows: (i) The oscillations are seen to decrease 
from above the quasiclassical H^f** 1 , due to the reentrant 
behavior of A^ Nc \ to be reduced considerably around 
B ~ 0.8ff c q 2 uasi , i.e. H*P S, /B ~ 1.25. However, they do 
not disappear completely in lower fields, (ii) This ex- 
tra attenuation is due to the broadening of the Landau 
levels caused by the pair potential, as in the case of the 
quadratic dispersion. Indeed, we have obtained quasipar- 
ticle spectra similar to those of Fig. [| (iii) The period of 
the oscillation remains unchanged above ~0.877^ asl , but 
some irregularity appears in lower fields. These features 
are in agreement with the results by Norman et alEB (iv) 
Little difference can be seen between the s- and d-wave 
attenuations. 



D. Summary of Two-Dimensional Calculations 

Let us summarize results and conclusions from our two- 
dimensional calculations, (i) Combining Figs. ||and|^, we 
are now able to attribute the origin of the extra dHvA 
oscillation damping unambiguously to the Landau-level 
broadening due to the scattering by the pair potential, 
(ii) As may be realized from Fig. o, presence of point 



nodes along the extremal orbit does not weaken -the at- 
tenuation, contrary to the statement by Miyake.Ej This 
fact suggests that the attenuation is determined by the 
average gap along the extremal orbit, (iii) The mean- 
field theory predicts that the dHvA oscillation comes to- 
gether with the oscillatory behaviors in H c2 and A^ N °', 
This will be so in three dimensionaljaodels where H C 2 (T) 
also shows an oscillatory behavior .C3 However, such sin- 
gular behaviors of H C 2 have never been identified defi- 
nitely in any materials displaying the dHvA oscillation, 
and reported H C 2 curves show more or less the smooth 
quasiclassical behavior. This discrepancy between the 
mean-field theory and the dHvA experiments remains 
a puzzle to be resolved in the future, (iv) The oscilla- 
tion attenuates considerably around B ~ 0.8iJ q 2 UaS \ i- e - 
H^2 asi /B^ 1.25, although we have set Huj h quasi/ fee T c ~l 

and Np 3> 1 at H^ 8 * 1 . Thus, the two dimensional mocU 
els fail to explain the experiment by Terashima et al.u 
which shows a persistent oscillation down to Q.2H C 2. In 
addition, the models cannot say anything about whether 
presence of a line node along the extremal orbit weakens 
the attenuation, (v) The approximation of retaining only 
A' ' works excellently for calculating M osc . Indeed, we 
have checked that the curves of M osc thereby obtained 
are almost indistinguishable from those of Fig. ||. (vi) 
The discrepancy mentioned in (iii) above suggests that 
we should rather use A^ ) obtained quasiclassically to re- 
produce the smooth behaviors of H C 2 in real materials. 
Figure m plots curves of M osc calculated using quasiclas- 
sical A™, i.e. the dotted lines of Fig. |5|. The oscillations 
are seen more regular than those of Fig. ^, but the am- 
plitudes attenuate almost similarly and are reduced con- 
siderably around H^^/B^ 1.25. We hence realize that 
using the quasiclassical A^ ^ suffices for the theory of the 
oscillation damping. This statement is especially true in 
the low-field region ~ 0.2£f^ asl where A^ ) approaches 
to the quasiclassical behavior, as may be realized from 
Fig. § 
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IV. THREE-DIMENSIONAL CALCULATIONS 

Having clarified basic features of the dHvA oscilla- 
tion for two-dimensional models as well as the mecha- 
nism of the extra oscillation damping, we proceed to con- 
sider three-dimensional models with various gap struc- 
tures which are more relevant to real materials. Our 
purposes in this section are summarized as follows: (i) 
to clarify the connection between the extra dHvA oscil- 
lation damping and the gap anisotropy by numerical cal- 
culations; (ii) to obtain an analytic formula for the extra 
oscillation damping; (iii) to estimate the gap magnitudes 
of various materials using the obtained analytic formula. 



A. Model 

The one-particle Hamiltonian (||) yields a spherical 
Fermi surface in the normal state. As for the pairing 
interaction, we consider three different models: 



v w(Ow(C) 

V 1 W(t)p-tW(?)&>6 



-K 2 ) 



(31) 



Here W(£) is a cut-off function given by Eq. ([23|), p = 
(sin 9 p cos ip p , sin 9 p sin ip p , cos 9 P ) specifies a point on the 
Fermi surface, and c = (sin 9 C , 0, cos 9 C ) denotes the direc- 
tion of the crystal c-axis, in the coordinate frame where 
B || z. Again the latter two models of Eq. (31 ) are beyond 
the original spherical interaction (jlTj), but they are con- 
venient for the above-mentioned purposes. In zero field, 
Eq. (gjj) yield 



A W(f ) p • c ia 3 a 2 



s-wave 



d x 2_ y 2-wave , (32) 



p z -wave 



which denote the isotropic s-wave state, a three- 
dimensional d 



x —y 



-wave state with four point nodes in 
the extremal orbit perpendicular to B, and the p-wave 
polar state with a line node perpendicular to c, respec- 
tively. The corresponding A^ c p ™' ) in Eq. (Q) can be writ- 
ten as 



t(N c tt, 



A( N ^W(£) sii 



8m2 + 8n 



■ l<7 



(33a) 



COS WnCOS 



5 ml +5 m -l 



where £= h 2 (N r /l 2 B +p 2 z )/2m c -e F , 8 p = tan" 1 (™^) , 
and A( Nc ~) is defined by 



( Vr 



B N / p , 

ti() ^ 2 

B jftpt 



Vi 



cos^p-cos 6 C *^pf ) 



sin t/pi sin 



$(^,l) + ^(^c-l) 



(33b) 



Here we have adopted a normalization for A^^ different 
from Eq. ( |20| ) so that this quantity acquires a direct cor- 
respondence to the maximum gap Ao in Eq. fl3^). The 
factors i in the second and third cases stem from cos 2(p p 
and cos tp p in Eq. (|33) , respectively. 

The coefficients A^) = A^) (B, T) in Eq. @ com- 
pletely specify the pair potential, as a lready mentioned. 
Based on the reasoning given in Sec. IHD(vi), we here 



adopt a quasiclassical A^ Nc ^ rather than the fully self- 
consistent one. Then the dominant A^ near H C 2 fol- 
lows the mean-field square-root behavior to an excellent 
approximation: 



A (0) =a{l-B/H c2 ) 1 



/2 



(34) 



See the dotted lines in Figs. and |^, for example. In 
addition, other components A^ c>0 ) can be neglected for 
the r elevant region B > 0.1H C 2, as pointed out in Sec. 
iHD| (v). We hence use the lowest-Landau-level approxi- 
mation of retaining only A^ ). The coefficient a — a(T) in 
Eq. (Q) is determined by requiring that the maximum 
of 



^JdR JdrA(r 1 ,r 2 )e- ipr / h 



(35) 



be equal to Aq(1 — B/H c2 ), where Aq(T) denotes the 
maximum gap obtained from the weak-coupling theory. 
This procedure yields 



'0.5, 



(36) 



for all the three cases of Eq. (|33|). Substituting Eq. (|3J) 
into Eq. ( |l3a|) with the choice fiw D ~ 10A (T = 0), the 
off-diagonal elements of Eq. ( |l5| ) are fixed completely. 

The above non-self-consistcnt procedure has another 
advantage that we can choose Hlob=h c2 an d £f in Eq. 
( [Hf ) independently. We have set 



hit) 



k B T c at T = , 



(37) 



in accordance with fiw^/feTc = 1 ~ 3 and Np ^> 1 for 
relevant materials (see Table H below). Also, we have 
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FIG. 8: (a) The oscillatory part of the magnetization M OBC in 
the vortex state (blue line) as compared with the normal-state 
one (sky-blue line) for the s-wave model of Eq. (^) at T = 0. 
(b) The corresponding Dingle plot (points with error bars) 
as compared with various theoretical predictions; see text for 
details. 



chosen s-p in such a way that there are about 50 Lan- 
dau levels below £f for the extremal orbit at H C 2- Now, 
the matrix elements of Eq. ( |T5| ) are specified completely. 
Hexagonal, square, and hexagonal lattices are assumed 
for the three cases of Eq. (p3), respectively. 



Numerical Procedures 



The wavevector p z of 0<p z < 1.2pF is discretized into 
~ 1000 points with an equal interval. For each of them, 
we have diagonalize d Eq. ( |l5| ) with the same procedure 



as described in Sec. |III B 
stituted into Eq 

Eq. 



The obtained results are sub- 
to calculate the magnetization by 
All the calculations are performed at T = 0. 



C. Numerical Results 

We first focus on the dHvA oscillation of the s-wave 
model in Eq. ( |32| ) to clarify its basic features. Using our 
numerical data, we also test the applicability of various 
theoretical formulas presented so far. 

Figure ||(a) presents oscillation of the s-wave magneti- 
zation (blue line) as compared with the normal-state one 
(sky-blue line). With ft^u^ = ^bT c , the oscillation is ob- 
served to persist down to a fairly low field of H C 2 /£? < 1 .8, 
i.e., B > 0.55H C 2, which is smaller than 0.8H C 2 around 
which hwB becomes equal to the spatial average of the 
energy gap, Eq. (|35|). This is partly because the gap is 
smaller within the extremal orbit, as shown quasiclas- 
sically by Brandt et alx3 Indeed, Fig. ^ calculated at 
£> = 0.968-ff C 2 demonstrates that the dispersion for p z = Q 
is smaller than that for p z = 0.9pF- This tendency re- 
mains in the high- field region of B>0.5H C 2- 

It has become conventional to express this extra atten- 
uation in the vortex state by introducing an additional 
factor R s for the dHvA oscillation amplitude: 



exp 



— exp 



27T 2 fc B T A 



(38) 




FIG. 9: Quasiparticle dispersion in the magnetic Brillouin 
zone for the s-wave model at B = 0.968i? C 2- (a) p z =0; (b) 
p z =0.9pF. 



where the parameters r s and Ta are directly connected 
with the extra Landau-level broadening r s in the vortex 
state as r s = h/2r s = 7rfc B XA- The points with error 
bars in Fig. ||(b) shows lni? s as a function of 1/B, i.e. 
the Dingle plot, obtained by numerical differentiation. 
This extra damping at high fields shows the behavior 
oc 1—B/H C 2 in the logarithmic scale, but irregularity sets 
in around 0.55i/ C 2 where the oscillation disappears. We 
attribute this irregularity to the effect of the bound-state 
formation in the core region. 

The lines in Fig. ||(b) are the predictions from various 
theoretical formulas. Maki's formulalla reproduces the 
correct functional behavior oc 1 — B/H C 2 at high fields, 
but the prefactor is seen too large. The NMA formula,E3 
deduced from the two-dimensional self-consistent numer- 
ical results with Np ~ 10 at H C 2 , predicts a more rapid 
attenuation incompatible with our numerical data. One 
reason for this discrepancy may originate from the fact 
that their numerical data with Np ~ 10 at H C 2 are not 
appropriate for obtaining an analytic formula by fitting. 
Another may be attributed to the difference in dimen- 
sions. Indeed, the dHvA oscillation in three dimensions 
differs from that in two dimensions on the point that 
some finite region 5p z around the extremal orbit is rel- 
evant. Most of the Landau levels in the region do not 
satisfy the particle- hole symmetry with respect to ef, so 
that the effect of the pair potential becomes smaller than 
that in twp-dimensions. Another theory by Dukan and 
Tesanovic,E-3 which would predict R s = in the clean limit 
of T = 0, is also inconsistent with the data. 

The red line in Fig. ||(b) is due to our formula for the 
extra Dingle temperature: 



fc B T A = 0.5f(|A p | 



m h c 1 — B/H c2 
ireh B 



(39) 



which is derived in Appendix |C| based on the second-order 
perturbation with respect to the pair potential. Here 
(| A p | 2 ) C o denotes the average gap along the extremal or- 
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FIG. 10: (a) The oscillatory part of the magnetization M OS c 
in the vortex state (blue line) as compared with the normal- 
state one (sky-blue line) for the d-wave model of Eq. (^) at 
T = 0. (b) The corresponding Dingle plot (points with error 
bars) as compared with the theoretical prediction (p9|). 



bit at B — 0, and is the band mass. The numerical 
constant 0.5 stems from Eq. (|36|), and f is a dimension- 
less quantity characterizing the Landau-level broadening 
due to the pair potential. This unknown parameter T is 
determined by a best fit to the s-wave numerical data, i.e. 
the points with error bars in Fig. |](b). This procedure 
yields 

f = 0.125. 

We observe in Fig. ||(b) that Eq. j39|), which predicts 
the dependence oc 1 — B/ H C 2 for In R s , agrees with the 
numerical results. This formula will be seen below to 
reproduce other numerical data excellently without any 
adjustable parameters. ■— ■ 

A difference of Eq. ( ^ ) from Maki's formulate lies in 
the prefactor where the Fermi velocity vp is absent. In- 
deed, a dimensional analysis on the second-order pertur- 
bation tells us that the Landau-level broadening in the 
vortex state should be of order \A^°\B)\ 2 /Hlub, where 



i p| 2 )eo' 



,(1—B/H C 2) is essentially the aver- 
age gap along the extremal orbit. This leads to Eq. j39| ) 
except for the numerical constant. 

We now turn our attention to see how the presence 
of point nodes affect the dHvA oscillation. Figure |Io|(a) 
shows the oscillation of the (i-wave magnetization (blue 
line) as compared with the normal-state one (sky-blue 
line). Although the d- wave gap in Eq. (^) has four 
point nodes on the Fermi surface along the extremal or- 
bit, the damping is seen strong and not much different 
from the s-wave case. From this fact, we may conclude 
that it is the average gap along the extremal orbit which 
is relevant for the extra dHvA oscillation damping. Fig- 
ure ^(b) presents the corresponding Dingle plot (points 
with error bars), which is compared with the prediction 
of Eq. (|39|). The formula with the average gap (|A p | 2 ) eo 
reproduces the numerical result for H C 2/B < 1.8 excel- 
lently without adjustable parameters, thereby providing 
a strong support for the above statement. 

This d-wave result is in disagreement with Miyake's 
theory that point nodes-jn the extremal orbit should 
weaken the attenuation.113 Indeed, Miyake's theory is 
based on a semiclassical quantization for the expression 
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FIG. 11: Left figures: the oscillatory part of the magnetiza- 
tion Most; in the vortex state (blue lines) as compared with 
the normal-state one (sky-blue lines) for the p-wave model of 
Eq. ([32]) at T = 0. The crystal c axis, which is perpendicular 
to the nodal plane, is tilted from the field B by 8 C = (top), 
9c~n/6 (second), and Oc — n/A (bottom). Right figures: the 
corresponding Dingle plots (points with error bars) compared 
with Eq. (El). 



of the electron number N e at B = 0. Neither his starting 
point N e (B = 0) nor the use of the semiclassical quantiza- 
tion may be justified for describing the dHvA oscillation 
observed mainly near H C 2- 

We finally consider the p-wave model with a line node 
in Eq. (|3|) to double-check the applicability of Eq. (|3S|). 
The left figures in Fig. [ll] display the dHvA oscillation 
for the line-node model, where the crystal c-axis is tilted 
from the magnetic- field direction by 9 C — (top) , 6* c = 7r/6 
(second), and 9 C = 7r/4 (bottom). The damping is seen 
weakest in the top figure where the gap vanishes along the 
extremal orbit, but increases gradually as finite gap opens 
along the orbit for C = — > ir/4. These results indicate 
conclusively that the average gap along the extremal or- 
bit is relevant for the extra dHvA attenuation. However, 
the non-zero extra damping in the top figure implies that 
not only the extremal orbit alone but some finite region 
around it contributes to the extra damping. Theoreti- 
cally, this corresponds to the fact that we have to perform 
the Fresnel integral Jl^expf— i(\/2iTNF p z /pp) 2 ] dp z for 
obtaining the LK formula in the normal state. Our data 
show that this off-extremal-orbit contribution cannot be 
neglected in the case where the gap vanishes exactly 
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TABLE I: Parameters characterizing three superconductors displaying the dHvA oscillation in the vortex state, together with 
values of the average gap along the extremal orbit \J (|A p | 2 ) co estimated by Eq. (J3S|) - Here the symbol a and 7 are band 



indices. These values of r/i(jAp| 2 } co are to be compared-with the band- and/or angle-averaged quantity A(0) extracted 
specific-heat experiment,Ej a far-infrared measurement ,lj a tunneling experiment,E3 or a Raman-scattering experiment] 


1 


3m a 


Compound T c Qs 


) H c2 (T=0) (T)_, m b /m„ ruv Hc2r ^) humJl 


~T C V<|Ap| 2 )co (meV) A(0)j 






NbSe 2 7.2t| 
Nb 3 Sn 18.3 
YNi 2 B 2 C 14.3 
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and completely at the extremal orbit. However, this off- 
extremal contribution is expected to become less impor- 
tant where finite gap is present along the extremal orbit. 
The right figures in Fig. [ll] show the corresponding Din- 
gle plot (points with error bars), which is compared with 
the prediction of Eq. (39). Except for the weak damping 
of 9 C = due to the off-extremal-orbit contribution, the 
formula is observed to reproduce the numerical results 
excellently. 



ESTIMATION OF ENERGY GAP 



nodes presented by Izawa et aZ.EzI based on a thermal- 
conductivity measurement. 



Our calculations in Sec. IV C have clarified that (i) the 
gap anisotropy can be detected by measuring the extra 
dHvA oscillation damping in the vortex state, and (ii) 
Eq. d39| ) is particularly useful for this purpose. Using 
the formula, we finally provide quantitative estimations 
of the average gap along the extremal orbit for several 
materials displaying the dHvA oscillation in the vortex 
state. Table | summarizes parameters describing three 
relevant materials. These materials commonly have fairly 
high T c 's, and the ratio Hluh^/^bTc ranges from 1 to 3. 
These features seem to be basic conditions for observ- 
ing the dHvA oscillation in the vortex state. The values 
for y/ (|Ap| 2 ) co are obtained by applying Eq. (^) to the 
observed dHvA attenuation in the vortex state. In do- 
ing so, we have adopted as m& in Eq. ([59]) the values 
from dHvA experiments rather than those from band, 
calculations, as indicated by the theory of Luttinger.E2l 
For comparison, we have also listed tha-yalues A(0) es- 
timated by a_specific-heat experiment ,c3La far- infrared 
measurement llH a turuacling experiment ,lj or a Raman- 
scattering experiment^ Thus, A(0) is expected to repre- 
sent band- and/or angle-averaged energy gap. As seen in 
Table |, the two quantities coincide excellently for NbSe2 
and NbaSn, indicating uniformly opened gap in these ma- 
terials. On the other hand, y/ (|A p | 2 ) co = 1.5 for the a 
band of YNI2B2C is considerably-smaller than A(0) = 2.5 
from a specific- heat experiment E3 This fact implies that 
YNi2B2C have large band- and/pr angle-dependent gap 
anisotropy. Indeed, Bintley et al!tB have recently carried 
out a detailed dHvA experiment on this material, rotat- 
ing the field direction and observing the extra attenua- 
tion. They have reported a large angle dependence of the 
attenuation magnitude. They have also pointed out that 
their result is in agreement with the model with point 



VI. 



SUMMARY 



We have carried out the first three-dimensional numer- 
ical calculations on the dHvA oscillation in the vortex 
state for various gap structures. We have thereby clari- 
fied the relation between gap anisotropy and persistence 
of the oscillation. We have also derived an analytic for- 
mula for the extra dHvA attenuation in the vortex state. 

Our main results are given by Figs. |-]ll] and Eq. ©. 
Those figures indicate clearly that the extra dHvA at- 
tenuation in the vortex state is directly connected with 
the average gap along the extremal orbit at B = 0. The 
derived formula (^) have been shown to reproduce the 
numerical results excellently. Our theory attributes the 
origin of the extra dHvA damping to the Landau-level 
broadening caused by the pair potential. Hence the peri- 
odicity of the vortex lattice assumed here is almost irrel- 
evant, and the theory is applicable also to the cases with 
irregularity such as a random array of vortices. Using Eq. 
(|39|), we have estimated average gap amplitudes along 
the extremal orbit for NbSe 2 , Nb 3 Sn, and YNi 2 B 2 C. 
The results indicate presence of large gap anisotropy in 
YM2B2C. 

Thus, we have shown explicitly that the dHvA effect in 
the vortex state can be a powerful tool to probe the aver- 
age gap along the extremal orbit. Our results imply that, 
by rotating the field direction and observing the attenua- 
tion amplitude, we can obtain unique information on the 
band- and/or angle-dependent gap structure. Such an 
experiment iias recently been performed on UPd2A}3 by 
Inada et alHi and on YM2B2C by Bintley et aZ.,L3 and 
the latter group indeed has detected large gap anisotropy 
in the ab plane. Equation ( |39| ) will be useful in similar ex- 
periments for estimating band- and/or angle-dependent 
gap amplitudes. 
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APPENDIX A: THERMODYNAMIC POTENTIAL Eq. (|T|) is defined by 



The Luttinger-Ward thermodynamic potential corre- 
sponding to Eq. (P is given by£3 



^Trln 



H-ie n A 
At -H*-ie 



e fe„0+ o 

e" l£ " + 



^TrA^, (Al) 



where we have adopted a compact notation of using 
x = ra with A air72 (ri, r 2 ) — » A(xi, x%), etc., and Tr also 
implies both integration and summation over r and a, 
respectively. The quantity e n /H denotes the Matsubara 
frequency, and 0+ is an infinitesimal positive const ant . 
Now, Eq. ([|) tells us .that the first matrix in Eq. ( Al) 
can be diagonalized asE3 



H(x,x') A(x,x') 
tf(x,x') -H*(x,x') 



E 



u*{x) -v s {x) 
v* s (x) -u s (x) 



E s 
-E A 



u s (x') v s (x') 
-v*(x') -u*(x') 



(A2) 



Substituting Eq. (A2) into Eq. flAl| ), the first term on 
the right-hand side becomes 



fc R T 



E / da; {[kMI 2 ^ 0+ + |^(x)| 2 e- z "°+] 



x \n(E s -z n ) 
+ [\v.(x)\ 2 e*»°+ + \u s (x)\ 2 e-^°+]\n(-E s -z n )}, 

(A3) 

with z n = ie n . The summation over n are then trans- 
formed with a standard techniqueEH into a contour in- 
tegral just above and below the real axis, using f(z) = 
(e z / l " BT + l) _1 and /(— z) for the terms with e z ™°+ and 
e -z„o + ^ reS p ec tively. Considering the poles inside the two 
contours and using /[|u s (x)| 2 + |?; s (a;)| 2 ]da; = 1, we obtain 
Eq. (0). 



APPENDIX B: BASIS FUNCTIONS AND 
OVERLAP INTEGRALS 

We here present explicit expressions for the quantities 
appearing in Eqs. ( |l2] ) and (|lg|); see Ref. |34| for their 
detailed derivations. It should be noted that we here 
adopt the symmetric gauge (||) which is more convenient 
than the Landau gauge used in Ref. [54| Hence there is 
an extra factor due to the gauge transformation in every 
expression of the basis functions, such as e~ lx vl 2l B m Eq. 
(|l|) below. 

The basis function ipNka (N = 0, 1, 2, • • • ; a = 1, 2) in 



AA f /2 

n=-N f /2+l 

Y -ixy/2l^~(x-l^k y ~na lx ) 2 /2lg+i(a-l)nTT 



x, i ^/t.^-^r™ 1 * ), (bi) 



Ib 



where S = n^ATf, and Hn(x) = e x2 e x2 is the 

Hermite polynomial.il The basis function ip^ a for the 



iVq 



center-of-mass coordinates is obtained from Eq. (Bl) by 
putting k^q, a = l, and Ib^I c = Ib/V^ as 

M/2 

ijjff (r)= V] e i[lv(y+^/2)+na la ,(y+llq x -nai v /2)/ll] 
n= _jV f /2+l 

xc~ 



^*}* /l ^- H N ( X -^- nai * ). (B2) 



2 n N\^kS 



Finally, the basis function ipNm for the relative coordi- 
nates, which is conveniently chosen as an eigenstate of 
the orbital angular momentum operator l z , is given by 



2<Kl I \(N+m)\ 



(B3) 



where (= (x + iy)/\/2Z r with Z r = \/2fe, and L^\x) = 



j^e, x x~ m (j^) N e~ x x N+m is the generalized Laguerre 
polynomial satisfying N + m > O.Ea 

We next provide expressions of the overlap integrals in 
Eq. ([l3|). The first one, which was obtaind by Rajagopal 
and RyanJ^I is given by 

(N c N t \N!N 2 ) 



i\Ti!J\r 2 !iV c !JV r ! 



min(AT 1 ,JV ) 



E 



r_l)AT 2 -JV c +n 



n=rrLa.x(0,N c — N 2 ) 



n\ (JVi-n)! (iV c -n)! (7V 2 -7V c +n)! 



(B4) 



which satisfies 

{N c N x \NxN 2 ) = (-^{JVcJVrlJVaJVi) , (B5a) 



^ (NMlNt^i^NzlNcNr) = 6 N , Nc 5 N , Nt . (B5b) 
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The asymptotic expression of Eq. ( B4 ) for iV c <C N± 
is given bya 



(JV c JV r |JViJV 2 ) » S Nl+N2 , Nc+Nr (-if 2 
vp -* 2 /2 H Nc (x) 



-,1/4 



7r(7V c +7V r ) 



(B6) 



with x=(N 1 -N 2 )/ y / 2(N c +N 1 .). This expression forms 
the basis for quasiclassical approximations. 
The second overlap integral is given by 

(q|m+iV r ) = V^l t (-1)™+"' [^Uq(°)]* ' ( B7 ) 

where ip^+N q is another basis function of the relative 
coordinates defined by 

Ms/ 4 

^W( r )= V" e i[g H (a+i?9x/4)/2+2na lx (y+i r 2 9x /2-na l!; )/i^ 
n=-M/4+l 



xe 



-iiry/2Z 2 -(:!;-2 2 ij s /2-2nai x ) 2 /22 2 



/ 2a lx /l T fx-l^q y /2~2na lx 



(B8) 



APPENDIX C: ANALYTIC FORMULA FOR 
DHVA OSCILLATION DAMPING 

In order to derive an analytic expression for the dHvA 
oscillation damping in the vortex states, we start from 
the thermodynamic potential of Eq. (||). The last term 
— 7;TrA<I> may be expressed solely with respect to the 
pair potential, so that they can be neglected in the 
present model to consider the oscillatory part. Since we 
are interested in the extra damping in the vortex state, 
we adopt as E s and v s the expressions from the second- 
order perturbation with respect to A. They are obtained 
as 

ENka Pz <j = I^JVp^l+^kp^sign^jvp^o-) , (CI) 



JWNk aPz cr(r)\ 2 dr = 6 , (-^jv P , CT )+?7^ k p 2 Sign(^Arp 2(T ) , 

(C2) 

where 6 is the step function, and is defined by 

using Eq. (13a) as 



(") = \ " 
7 ?JVkp z — / , 



|A 



(kp»)|2 
NN' I 



N 1 



(£iVp*+£iV'-pJ r 



(C3) 



The first terms on the right-hand side of Eqs. (CI) and 
(|C2|) are just the normal-state results. The second terms, 



on the other hand, denote the finite quasiparticle disper- 
sion in the magnetic Brillouin zone and the smearing of 
the Fermi surface, respectively, due to the scattering by 
the growing pair potential. It is useful to express 

in terms of M Q \B) and the cyclotron energy Hu>b of the 
extremal orbit as 



r ?A f kp i 



The quantity fj^„ thus defined is dimensionless, and 
we realize that the main B dependence in Eq. (C4) lies 
in the prefactor |A'°'(B)| 2 /(fiuB)". The explicit expres- 
sion of is obtained by using Eqs. (13a) and (|3^). 
Considering the case C = Q, for simplicity, it is given by 

~(n) _ AA f 2 x ^ \(NN'\0N+N')\ 2 (N+N' + m\2\<L-cti 



|AW(B)l 2 -(„) 
{faj B ) n VNkp * 



(C4) 



4 ^ 

N'mm 1 



[N+N'-2{N F + 6)Y 

SmoSm'o 



x(2k-q|iV + 7V'- 



5 m ,±2&rn' ,±2 Sm 4 p 

S m oS m >o cos 2 # p 



(C5) 



where the quantity 8 = 5(B,p z ) (\S\ < 1/2) specifies the 
location of e-p between the two closest Landau levels, and 
the overlap integrals are defined by Eqs. ( |B4|) and (|B7|). 
The corresponding normalized density of states: 



f ka 

will play a central role in the following. 
Substituting Eqs. (JClT) and flC^) into Eq. 



(2) 

that the terms containing rj^^ may be neglected due 
to the cancellation between the particle and hole contri- 
butions. The remaining, term can be transformed with 
the standard procedure.tl We thereby obtain, for the first 
harmonic of f2/V, the expression: 



(C6) 



we find 



B-l \ ^ 



pOO pOO /' CO 

' dN cos(2tt7V) / dp z / dfj 

-1/2 J-oc J-oo 



xD^(r?)ln[l+e-^-+^ A(0) ( B )l 2 /^ B )/ fcBT ] . (C7) 

The function D$ (fj) depends on (N,p z ), but may be 

replaced by a representative one-.D^ (fj) to be placed 
outside the N and p z integrals,^! where the recovered 
index i specifies the s-, d-, or p-wave case of Eq. jCq). 
It may also be acceptable to use a Lorenzian for it: 

Df ) (j7) = f^/7r(j? 2 + f|).@ We thereby obtain an expres- 
sion for the magnetization which carries an extra damp- 
ing factor: 



/oo 
Wfflexp -2mfj\A^(B)\ 2 /(Hiu B ) 2 dfj 
- OO 

= exp[-27rf,|A<°>(B)| a /(&«>B) 2 ] . (C8) 
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Thus, the superconductivity gives rise to an extra Dingle 
temperature of k B T A = f e \A^°\B)\ 2 /nhu} B , or equiva- 
lently, the extra scattering rate of t s _1 = 27t£:bTa/7i. 

Equation (CS) has an advantage that one can trace 
the origin of the extra dHvA damping definitely to the 
growing pair potential, which brings finite quasipa rtic le 
dispersion in the magnetic Brillouin zone as Eq. ( |C5| ), 
and the corresponding Landau-level broadening as Eq. 
(C6). Moreover, Eq. (C5) reveals that this broaden- 
ing near H c2 is closely connected with the zero-field gap 
structure given by Eq. (p2h. 



There seems to be no analytic way to estimate f s , so 
we fix it through the best fit to the numerical data of 
Fig. ||(b). Using Eq. (|||) with ci 2 = 0.5Aq, the procedure 
yields T s — 0.125, as mentioned before. It is also clear 
both from Figs. |o| and |l| and from Eq. ( |C5| ) that the 
average gap around the extremal orbit is relevant for the 
extra attenuation. We hence put a 2 T^ = 0.5(| A p | 2 ) eo r s . 
We thereby obtain Eq. (|3^), which yields excellent fits to 
the d- and p-wave numerical data without any adjustable 
parameters, as seen in Figs. |l^ and [ll[ 
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The p z inte gral in Eq. (G7) passes through regions where 
S~0 in Eq. (C5), i.e., the second-order perturbation is not 



justified. However, the main contribution to the integral 
cert ainly comes from the regions where we can use Eq. 
(CI). Notice that 6 = —1/8 at the maximal orbit when the 
magnetization takes a local maximum. 



It is found numerically that the variance of Eq. (C6) de- 
pends little on the values of Nf and N (~ Nf), as expected. 
It is ~0.2 for the s-wave pairing with 5 =—1/8. 



